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Abstract 

Image analysis problems, posed mathematically as variational principles or as partial differential 
equations, are amenable to numerical solution by relaxation algorithms that are local, iterative, 
and often parallel. Although they are well suited structurally for implementation on massively 
parallel, locally-interconnected computational architectures, such distributed algorithms are seriously 
handicapped by an inherent inefficiency at propagating constraints between widely separated 
processing elements. Hence, they converge extremely slowly when confronted by the large 
representations necessary for low-level vision. Application of multigrid methods can overcome 
this drawback, as we established in previous work on 3-D surface reconstruction. In this paper, we 
develop efficient multircsolution iterative algorithms for computing lightness, shape-from-shading, 
and optical flow, and we evaluate the performance of these algorithms on synthetic images. The 
multigrid methodology that we describe is broadly applicable in low-level vision. Notably, it is an 
appealing strategy to use in conjunction with rcgularization analysis for the efficient solution of a 
wide range of ill-posed visual reconstruction problems. 
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1. Introduction 

Variational principles and partial differential equations have played a significant role in 
the mathematical formulation of low-level visual information processing problems (representative 
examples include [Horn, 1974, 1975; Ulltnan, 1979; Horn & Schunck, 1981; Ikeuchi & Horn; 
1981; Narayanan el al, 1982; Bajcsy & Broit, 1982; Hummel & Zucker, 1983; Grimson, 1983; 
Terzopoulos, 1982, 1983; Nagel, 1983; Hildreth, 1984; Brady & Yuille, 1984]). An attractive feature 
of variational and differential formulations (once discretized) is the possibility of computing the 
desired solutions by a popular class of numerical relaxation algorithms. These iterative algorithms 
require only local computations which can usually be performed in parallel by many locally 
communicating processors distributed in computational networks or grids. 

Local, parallel algorithms are appealing in the context of low-level vision [Rosenfeld el 
al, 1976; Ullman, 1979; Ballard el al, 1983]. At a certain level of abstraction they do not 
appear incompatible with the apparent structure of advanced biological vision systems. Moreover, 
they are ideally suited to implementation on massively parallel computers with numerous simple, 
locally interconnected processing elements. Such potentially powerful architectures will certainly 
proliferate, pending imminent advances in VLSI technology [Batcher, 1980; Hillis, 1981]. 

The desired solutions to many visual problems appear to possess certain global properties 
/*""n (consistency, smoothness, minimal energy, etc.), which are expressed formally by the variational 

principle or associated partial differential equation formulations. 1 Given only local communication 
capabilities among processing elements, however, global properties can only be satisfied indirectly, 
typically by iteratively propagating visual constraints across the grid network. Indirect propagation 
can result in substantial computational inefficiency, since the computational grids necessary for low- 
level vision applications tend to be extremely large. Convergence of the iterative process is often 
so slow as to nearly neutralize the computational power offered by massive parallelism. Indeed, 
for fine discretizations on large grids, excruciatingly slow convergence rates have been observed in 
iterative algorithms for computing lightness [Blake, 1984; see also Horn, 1974], shape-from-shading 
[Ikeuchi & Horn, 1981; Smith 1982], optical flow [Horn & Schunck, 1981; Nagel, 1983], 3-D 
surfaces [Grimson, 1983; Terzopoulos, 1982, 1983], and other visual reconstruction problems. 

Since spaUal locality of computation is dependent on spatial resolution, local (e.g., nearest 

neighbor) computations on a coarse grid over a given region are analogous to more global 

computations on a fine grid over the same region. This suggests the possibility of counteracting 

die sluggishness of global interactions by deploying local iterative processes over a multiresolution 

hierarchy of grids. This is the basis of multigrid relaxation methods which are gaining popularity 

in applied numerical analysis [Hackbusch & Trottenberg, 1982]. The computational structure of 

1 Variational and differential formulations can be related through the Euler-Lagrange equations of the calculus 
^f**N of variations, given appropriate continuity and symmetry (or self adjoinuiess) conditions [Courant & Hilbert, 

1953]. 
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multigrid methods bears an interesting analogy to the multiresolution nature of spatial frequency 
channels in the human early visual system [Braddick, el al, 1978], The methods are also related 
to certain multiresolution image processing staictures that have been proposed, notably pyramids 
[Rosenfeld, 1984]. 

In earlier work, we developed an efficient surface reconstruction algorithm based on multigrid 
relaxation methods [Terzopoulos, 1982, 1983] and we suggested, as has Glazcr [1984], that multigrid 
methods are broadly applicable in low-level computer vision. After a brief overview of multigrid 
methodology, we apply it to three other vision problems: die well-known problems of computing 
lightness, shape-from-shading, and optical flow from images. We develop novel multiresolution 
algorithms for each problem. Our empirical results indicate that these algoridims offer order-of- 
magnitude gains in efficiency over their conventional single level counterparts. 

2. Multigrid Methodology 

Pioneering investigations into multigrid methodology include the work of [Fedorenko, 1961], 
[Bakhvalov, 1966], [Brandt, 1973, 1977], and [Nicolaides, 1977]. It has been applied to many 
boundary value problems (see [Brand, 1982] for an extensive bibliography) and tfiere has also been 
some development in the context of variational problems [Nicolaides, 1977; Biandt, 1980]. 

2.1. Multigrid Relaxation Methods 

Multigrid relaxation methods take advantage of multiple discretizations of a continuous 
problem over a range of resolution levels. The coarser levels trade off spatial resolution for 
direct communication paths over larger distances. Hence, they effectively accelerate the global 
propagation of information to amplify die overall efficiency of the iterative relaxation process. 

The inherent computational sluggishness of local iterative algorithms can be studied from a 
spatial frequency perspective. A local Fourier analysis of the error function (or, more conveniently, 
the dynamic residual function) from one iteration to the next shows that high-frequency components 
of die error — those components with wavelengths on die order of the grid spacing — arc short- 
lived, whereas low-frequency components persist through many iterations [Brandt, 1977]. Hence, 
common (L 2 or Loo) error norms decrease sharply during the first few iterations, so long as there are 
high-frequency components to be annihilated, but soon degenerate to a slow, asymptotic diminution 
when only low-frequency components remain (see Fig. 1). This suggests diat while relaxation is 
inefficient at completely annihilating die error function, it can be very efficient at smoothing it. 
From this point of view, the grid hierarchy enables die efficient smoothing properties of relaxation 
to be exploited over a wide range of spatial frequencies. 

Empirical studies of model problems (Poisson's equation in a rectangle) indicate that multigrid 
methods can converge in essentially order O(N) number of operations, where TV is the number of 



TERZOPOULOS 



MUI..TIGRID RELAXATION METHODS 




/-\ 



Figure 1. Asymptotic error reduction by relaxation. The mean square (dynamic residual) error is plotted 
as a function of the iteration number for a sequence of (Gauss-Seidel) relaxation iterations of a surface 
reconstruction algorithm. The curve exhibits a typical behavior of local iterative methods: Converger ce is 
rapid during the first few iterations, but quickly degenerates to slow asymptotic error reduction. 



nodes in the grid [Brandt, 1977]. This can be compared to typical complexities of 0(N 3 ) operations 
for the solution of model problems by standard (single level) relaxation. As a consequence, 
multigrid methods potentially offer dramatic increases in efficiency over standard relaxation methods 
in low-level vision applications, since N tends to be very large (order 10* to 10 6 , or more). For' 
comparative complexity analyses, the total computational expense of multigrid methods may be 
measured in convenient machine independent units. The basic work unit is defined as the ampunt 
of computation required to perform one iteration on the finest grid in the hierarchy. 

Our adaptation of multigrid methods to visual processing has a number of features: (i) 
multiple visual representations covering a range of spatial resolutions, (ii) local, iterative relaxation 
processes that propagate constraints within each representational level, (iii) local coarse-to-jfine 
prolongation processes that allow coarser representations to constrain finer ones, (iv) fine-to-coarse 
restriction processes that allow finer representations to constrain and improve die accuracy of coarser 
ones, and (iv) (recursive) coordination schemes that enable the hierarchy of representations and 
component processes to cooperate towards increasing efficiency. 

In multigrid metiiods, the intralevel processes usually are basic relaxation methods suci as 
Gauss-Seidel or Jacobi relaxation, the prolongation processes are local Lagrange (polynomial) 
interpolations, and die restriction processes are local averaging operations. The exact form of these 
operations is problem-dependent. 

2.2. Discretization 

Appropriate relaxation processes can be derived by local discretization of die continuous 
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problems. The finite element method [Strang & Fix, 1973], a general and powerful local 
discretization technique, can be applied directly to variational principle formulations of visual 
problems [Terzopoulos, 1982]. When die visual problem is posed as a partial differential equation, 
local discretization may be carried out using the finite difference method [Forsythe & Wasow, I960]. 

The basic idea behind the finite element method is diat a global approximation can result from 
interactions among many very simple local approximations. This is accomplished by tessellating 
the continuous domain into a large number of small subdomains or elements E whose dimensions 
depend on a fundamental size h. The approximation within elements depends on a small number of 
parameters — the values of the solution, and/or some of its derivatives, at a set of nodes associated 
with each element. The power of the method stems from the fact die local approximations can be 
based on low-order polynomials. This makes it relatively easy to express the continuous functional 
as a discrete summation over all the element contributions. If die variational principle is quadratic, 
the resuldng discrete problem takes the form of a large system of linear equations A h u h = f h , 
where u h is die vector of nodal variables. The finite element mediod can also be characterized as a 
systematic procedure for generating finite element approximating spaces whose local-support basis 
functions make A h sparse (i.e., most of its elements are zero). 

The finite difference method is applied differently. Typically a grid of nodes with spacings 
proportional to a parameter h is set up over the domain. The differential operator is then replaced 
/""""N by finite difference equations involving nodal variables at neighboring nodes. The collection 

of finite difference equations defines a discrete system which approximates the given differential 
equation. If the differential operator is linear (as are the Euler-Lagrange equations of quadratic 
variational principles) and a linear finite difference approximation is employed, the discrete system 
is again a linear system h h \x h = f 1 . Although die total number of nodes N is generally large, each 
finite difference equation involves only a few nodal variables. Therefore, the linear system is again 
sparse. 

While the finite difference method is generally easier to apply, die finite element mediod offers 
a much sounder convergence theory, as well as a flexibility that allows the spatially nonuniform 
discretization of domains having complicated shapes. Nonetheless, both discretization techniques 
yield large, sparse systems of linear equations in a wide range of visual applications. A great deal 
of effort in numerical analysis has been directed to the solution of such systems, which turn out to 
be especially well suited for solution by local, parallel, iterative methods, particularly (he relaxation 
methods that we have been discussing. 

2.3. Multigrid Structure and Coordination 
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Our spatially uniform discretizations of the continuous visual problems in this paper will 
yield uniform grids at each level of die multigrid hierarchy. Application of multigrid metiiods 
can be simplified substantially given a 2:1 decrease in grid resolution from any level to the next 
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Figure 2. Possible grid organization of a multiresolution algorithm. A small portion of three levels of the 2:1 
multigrid hierarchy is shown. Only nearest-neighbor int erprocessor connections ate included. 
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coarser level. Fortunately, this resolution ratio appears to be near optimal with regard to multigrid 
convergence rates [Brandt, 1977]. Fig. 2 illustrates a portion of three grids of a 2:1 multigrid 
hierarchy. In a serial implementation the central processor operates at each grid node sequentially, 
whereas in a fully parallel implementation, each node represents a scpkrate processing element 
within a distributed local-interconnect architecture (see Fig. 2). j 

The multiresolution visual algorithms to be described utilize simple injection I (=> ,_i for 



the fine- to-coarse restrictions, bilinear interpolation I ( _ 1=> ( for the cor 
and an adaptive multigrid coordination scheme which was employed su 
reconstruction algorithm (see [Terzopoulos, 1982, 1983] for details) 



rse-to-fine prolongation, 
:cessfully in our surface 
'he general coordination 



scheme first performs a sufficient number of relaxation iterations to solve the coarsest level discrete 
system A ;il u> = f* 1 to desired accuracy (procedure SOLVE), and then proceeds to the finest level 
/ = L according to 



/"""N 



/"""S 



TERZOPOULOS 




procedure FM6 

u' 1 ' - SOLVE (1, u\ 
for / <- 2 to Z, do 
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MG {I, V h ', f") 
end; 

applying the multigrid algorithm 

procedure MG (/, u, g) 

if 1 = 1 then u «- SOLVE (1, u, g) 
else 
begin 

for i <~ 1 to ni [while ...] do u +- RELAX (I, u, g); 
v«- Ij=>/-iii; 

d<- A^-'v+I^i-ifg-A^u); 

for i <- 1 to n 2 [while ...] do MG {l-l, v, d); 
u <- u + I/_i=*;(v - Ii=>(_i u) ; 

for i «- 1 to n 3 do [while ...] u «- RELAX (/, u, g) 
end; 

After ni relaxation iterations (procedure RELAX) have been performed at level /, MG performs a 
restriction to the next coarser level I — 1. It men calls itself recursively on the coarser level n 2 
times. Finally, it performs a prolongation from the coarser level back to level I, following up with 
n 3 more iterations on level I. The equations on the coarsest level I — 1 may be solved to desired 
accuracy with sufficiently many iterations (procedure SOLVE). One can readily show that when MG . 
is invoked on level X it calls RELAX a total of ^"'(nx + n 3 ) times on level I =^ 1 and it calls 
SOLVE n 2 x_1 times on level 1. In general, most of the relaxation iterations are performed on the 
coarser levels [Hemker, 1980]. 

The optional [while . . .] clauses denote conditions that may be checked during the 
computation and used to terminate some iterations. Dynamic conditions, typically convergence 
rates measured by error norms, are incorporated into adaptive coordination schemes, whereas fixed 
schemes are controlled only by the constants n x , n 2 , and n 3 [Brandt, 1977]. Although adaptive 
schemes tend to be more efficient in practice, fixed schemes lend themselves better to theoretical 
analysis and, moreover, they are easier to implement on distributed local-interconnect architectures 
due, in part, to the absence of error norms which require global computations. 



3. The Lightness Problem 

The lightness of a surface is the perceptual correlate of its reflectance. Irradiance at a point 
in the image is proportional to the product of the illuminance and reflectance at the corresponding 
point on the surface. The lightness problem is to compute lightness from image irradiance, without 
/""*\ any precise knowledge of either reflectance or illuminance. 
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3.1. Analysis 

The rctinex theory of lightness and color proposed by Land and McCann [1971] is based on the 
observation that illuminance and reflectance patterns differ in their spatial properties. Illuminance 
changes are usually gradual and, therefore, typically give rise to smooth illumination gradients, 
while reflectance changes tend to be sharp, since they often originate from abmpt pigmentation 
changes and surface occlusions. Horn [1974] proposed a two-dimensional generalization of the 
Land-McCann algorithm for computing lightness in Mondrian scenes, consisting of planar areas 
divided into subregions of uniform matte reflectance. 

Let R(x, y) be the reflectance of the surface at a point corresponding to the image point (x, y) 
and let S(x, y) be the illuminance at that point. The irradiance at the image point is given by 
E(x, y) = S(x, y) x R{x, y). Denoting the logarithms of the above functions as lowercase quantities, 
we have e(x, y) — a(x,y) + r(x,y). Applying the Laplacian operator A gives d(x, y) = Ae(x,y) = 
As(x, y)+Ar(x, y). In a Mondrian, illuminance is assumed to vary smoothly so that As(x, y) is finite 
everywhere, while Ar(x, y) exhibits pulse doublets at intensity edges separating neighboring regions. 
A thresholding operator T can be applied to discard the illuminance component: T[d(x, y)] = 
Ar(x,y) = f(x,y). Hence, the reflectance R is given by the inverse logarithm of die solution to 
Poisson's equation 

/**■% Ar(x, y) — f(x, y), in fi, 

where fi is the planar region covered by the image. 

Horn solved die above partial differential equation by convolution with the appropriate 
Green's function. We instead pursue a local, iterative solution based on the finite difference 
method. Suppose that fi is covered by a uniform square grid with spacing h. We can 
approximate Ar = r xx + r yy using the order h 2 approximations r xx = (r-* +liJ - - 2r£y + r^_ 1;j )/h 2 
and r yy = (r£y +1 - 2r£ y + r^J/Zi 2 to obtain a standard discrete version of Poisson's equation 
(if+i,y + i'i-ij + r?,y+i + ti,j-i - 4r^_,)//i 2 — f^-. This denotes a system of linear equations with 
sparse coefficient matrix. 

Rearranging, the Jacobi relaxation step is given by 

h (n+l) _ 1/fc (n) h (n) h (n) h (n) 2fh \ 

T i,3 — JV'i+l-i +T '-t,3 +r «',i+l +r t,i-l ~ nl i,j)> 

where the bracketed superscripts denote the iteration index. Jacobi relaxation is suited to parallel 
synchronous hardware, whereas the Gauss-Seidel relaxation step given by 

r h («+!)_ J /> W, r h (n+l), r h ("),,.* ("+ 1 ) u2fh \ 
r «,J — 4^«'+l. 3 +T i~i,3 +T i,3+l +l i,3-l ~ nl i,j) 

is more suitable on a serial computer and, moreover, requires less storage. 
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Figure 3. Synthesized Mondrian images. These images, input to the algorithm, contain patches of uniform 
reflectance and a left-to-right illumination gradient. The three smaller images are increasingly coarser sampled 
versions of the largest image which is 129 x 129 pixels, quantized to 256 irradiance levels. 

We note ia passing that Poisson's equation Ar = / is the Euler-Lagrange equation for the 
variational principle associated with a membrane problem. The solution can be characterized as 
the deflection v(x, y) — r(x, y) of a membrane subject to a load f(x, y), and it minimizes the 
potential energy functional £(v) = f f n \{v\ + v\) - fvdxdy [Courant & Hilbert, 1953]. Blake 
[1984] offers an alternative variational principle for lightness. Posing the lightness problem as 
a variational principle permits the direct application of the finite element discretization method, 
which for instance does not require a uniform discretization of fi. 

3.2. Results 



/ #- N 



A four level multiresolution lightness algorithm (with grid sizes 129 x 129, 65 x 65, 33 x 33, and 
17 x 17) was tested on a synthesized Mondrian scene consisting of patches of uniform reflectance, 
subjected to an illumination which increases quadratically from left to right. The original image, 
which is 129 x 129 pixels in size, and three coarser-sampled versions are shown in Fig. 3. All 
images are quantized to 256 irradiance levels. The grid function f*y, shown in Fig. 4, was 
computed by maintaining only the peaks in the Laplacian of r^-. Zero boundary conditions were 
provided around the edges of the images, and the computation was started from the zero initial 
approximation r^. '— o. 
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Figure 4. The grid function f?, y on each level. These functions were obtained by maintaining only the peaks 
hi the Laplacian of r£y at each level. 

- ' ' * ! 

Fig. 5 shows the reconstructed Mondrian which now lacks most of the illumination gradient. 
Reconstruction of the image from the functions shown in Fig. 4 required 33.97 work units. The 
total number of iterations performed on each level from coarsest to finest respectively is 142, 100, 
62, and 10. In comparison, a single- level lightness algorithm required about 500 work units to 
compute a solution of the same accuracy at the finest level in isolation. The single-level algorithm 
requires at least as many iterations for convergence as there are nodes across the surface, since 
information at a node propagates only to its nearest neighbors in one iteration. The multilevel 
lightness algorithm is much more efficient because it propagates information more effectively at the 
coarser scales. 
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4. The Shape-From-Shading Problem 

In general, image irradiance depends on surface geometry, scene illuminance, surface 
reflectance, and imaging geometry. The shape-from-shading problem is to recover the shape of 
surfaces from image irradiance. By assuming that illuminance, reflectance, and imaging geometry 
are constant and known, image irradiance can be related directly to surface orientation. 

4.1. Analysis 

Let u(x,y) be a surface patch with constant albedo defined over a bounded planar region fi. 
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Figure 5. The reconstmcted Mondrian. This is the solution computed after 33.97 work units by the four-level 
lightness algorithm. Most of the illumination gradient in Fig. 3 has been eliminated. 

The relationship between the surface orientation at a point (x,y) and the image irradiance there 
E(x, y) is denoted by R(p, q), where p = u x and q == u y are the first partial derivatives of the surface 
function at (x, y). The shape-from-shading problem can be posed as a nonlinear, first-order partial 
differential equation in two unknowns, called the image-irradiance equation: E(x, y) - R(p, q) — 
[Horn, 1975].. Surface orientation cannot be computed strictly locally because image irradiance 
provides a single measurement, while surface orientation has two independent components. The 
image irradiance equation provides only one explicit constraint on surface orientation. 

Ikeuchi and Horn [1981] proposed an additional surface smoothness constraint and the 
use of surface occluding contours as boundary conditions. Since the p-q parameterization of 
surface orientation becomes unbounded at occluding contours, however, surface orientation was 
reparameterized in terms of the (bounded) stereographic mapping: / = 2ap, g = 2aq, where 

a = 1/(1 + y/l + p2 + ? 2). 

These considerations are formalized by a variational principle involving die minimization of 
the functional 

£(f,9)^fJ n (rt + fl) + (9l + 9l)dxdy + ±JfjE(x, y )-R(f,g)} 2 dxdy. 

The first integral incorporates the surface smoothness constraint. The second is a least-squares term 
which coerces die solution into satisfying the image irradiance equation by treating the equation as 
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Figure 6. Lambertian sphere images. These synthetic images input to the algorithm show a Lambertian 
sphere distantly illuminated from the viewing direction. The three smaller images are increasingly coarser 
sampled versions of the largest image which is 129 x 129 pixels, quantized to 256 irradiance levels. 

a penalty constraint weighted by a factor X. Other variational formulations for shape-from-shading 
have been suggested, e.g., [Brooks & Horn, 1984]. 

The Euler-Lagrange equations are given by the system of coupled partial differential equations 

&f-\[E(x,y)-R(f,g)}R f =0, 
Ag-\{E(x,y)-R(f,g)}R g = 0. 

Discretizing these equations on a uniform grid with spacing h using standard finite difference 
approximations yields the Jacobi relaxation scheme 



f? 



(n+l) 



i.J 



jfffW^fc^ 



*[$,]« + X[2Su - R($J n \ g?,/ n) )][a,]JJ, 
8?./ n+l) - W } + X[JS?« - R($J n) , g?,/ B) )][iy&, 

where #[$,] = [f?_ 1>i + f? +1>y + ^-i + l? f y +1 ]/4 and #[ g y = [g*_ l|i + g? +lti + g? li _ 1 + g?, i+1 ]/4 
are local averages of Z' 1 and g H at node (i,y) (a factor of 1/4 has been absorbed into X), Rf = 
dR/df, and # 9 = dR/dg. On a sequential computer, we prefer to use the analogous Gauss-Scidel 
relaxation in our multilevel algorithm, due to its greater stability, faster convergence, and reduced 
memory requirements. Appropriate boundary conditions can be specified at occluding contours in 
the image. 
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Figure 7. Surface normals of the Lambertian sphere. The solution at the four resolutions that were obtained 
after 6.125 work units are shown. 

4.2. Results 
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A four level shape-from-shading algorithm (with grid sizes 129 x 129, 65 x 65, 33 X 33, and 
17 x 17) was tested on a synthetically-generated image of a Lambertian sphere distantly illuminated 
from the viewing direction by a point source. The original image, which is 129 X 129 pixels in size, 
and three coarser-sampled versions are shown in Fig. 6. All images are quantized to 256 irradiance 
levels. For the Lambertian surface, we employed the expression R{f,g) = max[0, cost], where 
cost == [16(/./ + g a g) + (4 - / 2 - <? 2 )(4 - f\ - g \))/[(A + / 2 + <7 2 )(4 + / 2 + g])] and where /. and 
g a are the light source direction components [Ikeuchi & Horn, 1981], and analogous expressions 
for its derivatives Rf and R g . The orientation of die surface was specified around the occluding 
contour of die sphere, and by treating die contour itself as a possible orientadon discontinuity, die 
grid functions / and g were allowed to make discontinuous transitions across it. Computation was 
started from the zero inidal approximation / — g = 0. 
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Figure 8. Surface representations of the Lambertian sphere. The depth representations on the left were 
generated by a four-level surface! reconstruction algorithm in 8.8 work units using the normal vectors in 
Fig. 7 as orientation constraints. On the right, the orientation constraints are depicted as "needles" on the 
reconstructed surfaces. Only the three coarsest levels are shown, since the finest resolution surface is too dense 
to render as a 3-D perspective plot. 
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The solution computed at the four levels after 6.125 work units are shown in Fig. 7. The 
total number of iterations performed on each level from coarsest to finest respectively is 32, 10, 4, 
and 4. In comparison, a single-level algorithm required close to 200 work units to obtain a solution 
of die same accuracy at the finest level in isolation. As in the case of the lightness problem, the 
single-level algorithm requires at least as many iterations for convergence as there are nodes across 
the surface, since information at a node propagates only to its nearest neighbors after each iteration. 
Convergence is somewhat faster, however, because shading information is available at every node 
inside the occluding contour to constrain surface shape according to die image irradiance equation. 
In any case, the multilevel shape-from-shading algorithm is again much more efficient because it 
enables information to propagate quickly at the coarser scales. 

To obtain a representation of the surface in depth, the surface normals in Fig. 7 were 
introduced as orientadon constraints to a four-level surface reconstruction algorithm with identical 
grid sizes [Terzopoulos, 1984a]. The normal vectors were first transformed from the f-g 
stereographic parameterization used in the shape-from-shading algoridim to die p-q gradient space 
parameterization used in the surface reconstruction algorithm using die formulas p — -4f/(f 2 + 
g 2 - 4) and q — -4g/(/ 2 + g 2 - 4). Nodes outside the occluding contour of the sphere were treated 
as depdi discontinuities. Fig. 8 shows die surfaces generated by the algorithm at the three coarsest 
resolutions. The reconstruction required an addidonal 8.8 work units. 

5. The Optical Flow Problem 

Optical flow is the distribution of apparent velocities of irradiance patterns in the dynamic 
image. The velocity field and its discontinuities can be an important source of information about 
the configurations and motions of visible surfaces. The optical flow problem is to compute a 
velocity field from a temporal series of images. 

5.1. Analysis 

Horn and Schunck [1981] suggested a technique for determining optical flow in the restricted 
case where the observed velocity of image irradiance patterns can be attributed directly to small 
interframe motions of surfaces in the scene. Under diese circumstances, the change in image 
irradiance at a point (x, y) in die image plane at time t and the motion of the irradiance pattern 
can be related by the flow equation E x u + E y v + E t ~ 0, where E(x, y,t) is the image irradiance, 
and u — dx/dt and v — dy/dt are the optical flow component functions. _ h 

An additional constraint is needed to solve this linear equation for the two unknowns u and v. 
If opaque objects undergo rigid motion or deformation, most points have a velocity similar to diat 
of dieir neighbors, except where surfaces occlude one another. Observing that the velocity field 
varies smoodily almost everywhere, optical flow can be determined by finding the flow functions 
u(x, y) and v(x, y) which minimize the functional 
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£(u, v) = a 2 J J {u\ + u\) + [y\ + v\) dx dy + J J {E x u + E y v + E t ) 2 dx dy, 

where a is a constant. The first term is the smoothness constraint, while the second is a least-squares 
penalty expression which coerces the flow field into satisfying the flow equation. Related variational 
formulations of the optical flow problem have been suggested (e.g., [Nagel, 1983], [Cornelius and 
Kanade, 1983]). 

The Euler-Lagrange equations for the functional £ are given by [Horn and Schunck, 1981] 

E\u + E x E y v = a 2 Au - E x E t , 
E x E y u + E\v — a 2 Av - E y E t . 

Assuming a cubical network of nodes with spacing h, where i, j, and k index nodes along the x, y, 
and t axes respectively, we use the following finite difference formulas to discretize the differential 
operators: 

[ E z]i,j,k — 2h,( E i+h3,k ~ E i-l,j,k)> 
\ E v\i,j,k = 2h,( E i,i+Uk - E i,j-l,k)> 

\ E t)i,i,k = ^{ E i,j,k+1 - E i,],k), 
A*— ^(*[u?,y,J -<,»), 

where *[\# Jtk ] = i(u?_ Wfk + u? ii+1]i + u? +ljiiJt + ufo_ 1)fc ) and 

^[ v ?,i,fc] = 4( v *-w,* + v ?,i+i,* + v i+t,i,* + vfli-i,*)- 0ther approximations are possible, including 
those suggested by Horn and Schunck which, however, require over four times the computation per 
iteration to gain some improved attenuation of high frequency error. Given dynamic images over 
at least three frames, a symmetric central difference formula [E t }^ >k = ^x (■#£;,*+! ~ E i,i,k-i) 
would be preferable, provided it is stable. 

Substituting the above approximations into the Euler-Lagrange equations and solving for 
Uij,k and v£y fc yields the following Jacobi relaxation formulas 

where rf dtk = {[E x ]^ k f + ([E,]^)* + £« 2 and 

"i,i,k = [■&]?,,•,*#["&,*] + [^v]?,y,**[ v w,*] + l E t]i,j,k- The natural boundary conditions of the zero 
normal derivative are appropriate on the boundaries of surfaces. They can be enforced by copying 
values to boundary nodes from neighboring interior nodes. 
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Figure 9. Lambertian sphere images. These synthetic images input to the algorithm at four resolutions depict 
a uniformly expanding Lambertian sphere, distantly illuminated from the viewing direction. Frames for the 
first time instant are shown to the left of frames for the secon d time instant 

5.2. Results 

A four level optical flow algorithm (with grid sizes 129 x 129, 65 x 65, 33 x 33, and 17 x 17) 
was tested on a synthetically-generated image of a Lambertian sphere distantly illuminated from 
the viewing direction by a point source. The sphere expanded uniformly over two frames. The 
first frame, which is 129 x 129 pixels in size, and three coarser-sampled versions are shown in the 
left half of Fig. 9. The next frame, in which the sphere has expanded is shown in the right half 
of die figure. All images are quantized to 256 irradiance levels. The velocity field was specified 
around the occluding contour of die sphere, and by treating the contour as a possible flow field 
discontinuity, « and v were allowed to make discontinuous transitions across it. The computation 
was started from the zero initial approximation u — v = 0. 

The solution computed on the three coarsest levels after 4.938 work units are shown in Fig. 
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Figure 10. Velocity vectors for the expanding Lambertian sphere. The solution at the three coarsest resolutions 
that were obtained after 4.938 work units are shown (the finest-level solution is too dense to depict). 

10 as velocity vectors in ay-space. The total number of iterations performed on each level from 
coarsest to finest respectively is 40, 5, 4, and 3. In comparison, a single-level algorithm required 
37 work units to obtain a solution of the same accuracy at the finest level in isolation. Again, 
the multilevel algorithm is more efficient because it propagates information quickly at the coarser 
scales. Glazcr [1984] also reports improvements consistent with ours with regard to the convergence 
rate of a multilevel optical flow algorithm relative to a single level algorithm. He employed the 
Horn-Schunck relaxation formulas for his implementation. 
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6. Multigrid Methods, Regularization, and Stochastic Relaxation 

A primary purpose of low-leyel visual processing is to reconstaict relevant physical charac- 
teristics of 3-D scenes from their images. We have considered in this paper three different visual 
reconstruction problems — the computation of lightness from an image (a 2-D, static reconstruction 
problem), shape-from-shading (a 3-D, static problem), and optical flow (a 2-D, dynamic problem). 
It was possible to apply multigrid mediods because each of these problems was formulated as a 
variational principle or associated partial differential equation. 

As inverse mathematical problems, visual reconstruction problems tend to be mathematically 
ill-posed, in that existence, uniqueness, and stability of their solutions cannot be guaranteed a 
priori [Poggio and Torre, 1984]. Among the systematic techniques diat have been developed to 
tackle ill-posed problems is the method of regularization [Tikhonov and Arsenin, 1977]. Through 
regularization analysis, ill-posed visual problems can be restated as well-posed variational principles 
by restricting the possible solutions with appropriate stabilizing functionals. In general, the 
smoothness properties of stabilizers must be controlled near discontinuities [Terzopoulos, 1984b]. 
Interestingly, the same stabilizer was used to impose the smoothness constraint in both the shape- 
from-shading and optical flow problems. 

A major attraction of regularization analysis is that it leads systematically to variational 
principles which permit advantageous use of multigrid relaxation methods. As a visual algorithm 
design strategy, regularization analysis applied in conjunction with multigrid methodology promises 
to impact on a broader spectrum of visual reconstruction problems, including image reconstruction 
and discontinuity detection [Geman and Geman, 1984], stereopsis [Marr and Poggio, 1977], 
registration [Bajcsy & Broit, 1982], motion field interpolation [Hildreth, 1984], shape-from-contour 
[Brady & Yuille, 1983], and structure-from-motion [Ullman, 1979]. 

An issue of concern is that the regularization of visual reconstruction problems cannot always 
be expected to lead to convex variational principles having a unique absolute extremum, without 
relative extrema. Unfortunately, classical relaxation or gradient descent mediods are not directly 
applicable to nonconvex variational principles, since they often get trapped in relative extrema. 
Stochastic relaxation algorithms (such as simulated annealing) do not suffer this disadvantage 
[Kirkpatrick el at, 1983; Hinton & Sejnowski, 1983]. Nonetheless, since stochastic relaxation 
searches for absolute extrema with processors that are restricted to local interactions, it too suffers 
serious inefficiencies in propagating constraints. The inherently slow convergence rates are further 
aggravated by the nondeterministic nature of the local computations. Multigrid methods may 
ameliorate these problems by facilitating constraint propagation through the use of coarser scales. 

7. Conclusion 

Many important problems in low-level computer vision can be formulated as variational 
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principles or as partial differential equations. A particular source of such formulations is the 
regularization analysis of ill-posed visual reconstruction problems. Once discretized, variational and 
differential formulations are amenable to numerical solution by iterative relaxation methods, which 
readily map into massively parallel computer architectures. However, distributed local-support 
computations are inherently inefficient at propagating constraints over the large network or grid 
representations that are encountered in computer vision applications. 

In our previous work on surface reconstruction algorithms, we established that multiresolution 
relaxation techniques can overcome this inefficiency, without sacrificing the local-interconnect nature 
of the computations. This has been corroborated in the present paper by successfully applying 
multigrid methods to the well-known problems of computing lightness, shape-from-shading, and 
optical flow from images. The novel multiresolution algorithms that we designed in the context 
of each of these problems were shown to be substantially more efficient than the published single 
level versions. 

Beyond its effectiveness as a (local) convergence acceleration strategy, our adaptation of 
multigrid methodology also leads to iterative algorithms that compute mutually consistent visual 
representations over a range of spatial scales. Multiresolution representations appear to be crucial 
in interfacing low-level visual processing to subsequent tasks such as recognition, manipulation, 
and navigation. 
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